* (c) M. Ciccarelli. Rats version 6.2
*
** program to do forecast using previously estimated components.
** compare the performance with RW AR PHIL. Reproduce Table 8.
**
** the model is the usual stock-watson
*
** y_h(t+h) = c + a(L)*y(t) + d(L)*F(t)+ e(t)
*
** where notation is as in the paper.
** In particular:
*
** 1. y_h(t+h) = (400/h)log(P(t)/P(t-h))
*
** 2. y(t) = 400*log(P(t)/P(t-1))
*
*
*
** The model is compared withthree benchmarks
*
** a. AR
** b. Phillips curve
** c. RW
*
*
OPEN DATA forecast_DATA.xls    ;* SA data
CALENDAR 1960 1 4
ALL 2008:02
DATA(FORMAT=XLS,ORG=COLUMNS) 1960:01 2008:02
close data
comp lags1 = 1
comp seas = 1

if seas ==1
{
X11(GRADUATE) CPIQU2 / CPIQU2
X11(GRADUATE) CPIQAT / CPIQAT
X11(GRADUATE) CPIQLU / CPIQLU
X11(GRADUATE) CPIQIT / CPIQIT
X11(GRADUATE) CPIQPT / CPIQPT
X11(GRADUATE) CPIQNL / CPIQNL
X11(GRADUATE) CPIQSE / CPIQSE
X11(GRADUATE) CPIQBE / CPIQBE
X11(GRADUATE) CPIQDE / CPIQDE
X11(GRADUATE) CPIQDK / CPIQDK
X11(GRADUATE) CPIQES / CPIQES
X11(GRADUATE) CPIQFI / CPIQFI
X11(GRADUATE) CPIQuk / CPIQUK
X11(GRADUATE) CPIQFR / CPIQFR
X11(GRADUATE) CPIQGR / CPIQGR
X11(GRADUATE) CPIQIE / CPIQIE
X11(GRADUATE) CPIQUS / CPIQUS
X11(GRADUATE) CPIQJP / CPIQJP
X11(GRADUATE) CPIQAU / CPIQAU
X11(GRADUATE) CPIQCA / CPIQCA
X11(GRADUATE) CPIQNZ / CPIQNZ
X11(GRADUATE) CPIQNO / CPIQNO
X11(GRADUATE) CPIQCH / CPIQCH
}
end if

comp n = 23
comp nsteps = 8
comp lags = 0

dec vect[series] ip(n) m3(n)
dec vec[series] infla(n)
dec vec[series] infla1(n)
dec vec[series] infla_steps(n)
dec vec[series] fore_infl(n)
dec vec[series] fore_infl_phil(n)
dec vec[series] fore_infl_ar1(n)
dec vec[series] fore_infl_rw1(n)
dec rec[series] rmse_far_serie(n,3) rmse_naiv_serie(n,3) rmse_ar1_serie(n,3) rmse_phil_serie(n,3)

comp lag_bic = 2;*  max lag to search with bic criterion
comp lagopt = 1;* 0 = no optimal choice of lag
comp lag_fact1 = 4
comp lag_fact2 = 2
comp lag_fix = 2
COMP AVER = 0

dec vec[integer] opt_lag(n)
dec vec bic(lag_bic)

source(noecho) princomp.src
source(noecho) uforeerr.src

set dcp = log(comprice/comprice{1})*400

comp beg = 1960:3
comp tempend = 1985:4
comp tempend1 = 1994:4
comp end = 2008:2
comp fin = end-tempend
comp wind = tempend1-tempend
comp count_st = 0

dofor steps =   1  4 8

 comp count_st = count_st + 1

 disp steps
*
 do i = 1,n
  set infla(i) = 400.*log(([series] i)/(([series] i){1}))
  set infla_steps(i) = (400./steps)*log(([series] i)/(([series] i){steps}))
 end

 comp J=1;
 dofor i = M3QU2 M3QUS M3QCA	M3QUK	M3QJP	M3QDE	M3QFR	M3QIT	M3QAT	M3QBE	M3QDK	M3QFI	M3QGR	M3QIE	M3QDE M3QPT	M3QES	M3QSE	M3QNL	M3QAU	M3QNZ	M3QNO	M3QCH
  set m3(j) = log(i{0}/i{1})*400
  comp J=J+1;
 end dofor

 comp J=1;
 dofor i = IPQU2	IPQUS	IPQCA	IPQUK	IPQJP	IPQDE	IPQFR	IPQIT	IPQAT	IPQBE	IPQDK	IPQFI	IPQGR	IPQIE	IPQLU	IPQPT	IPQES	IPQSE	IPQNL	IPQAU	IPQNZ	IPQNO	IPQCH
  set ip(j) = log(i{0}/i{1})*400.0
  comp J=J+1;
 end dofor

* smpl beg end

 clear forint FORE_INFL FORE_INFL_RW1 FORE_INFL_AR1 FORE_INFL_PHIL

* clear forint

 comp inum = 0

 DO ti = steps,fin
*
  comp count = 1
   dofor jj = infla(1) infla(2) infla(3) infla(4) infla(5) infla(6) infla(7) infla(8) $
    infla(9) infla(10) infla(11) infla(12) infla(13) infla(14) infla(15) infla(16)$
    infla(17) infla(18) infla(19) infla(20) infla(21) infla(22) infla(23)
    stat(noprint) jj beg tempend+ti-steps
    set infla1(count) = (([series] jj)-%mean)/sqrt(%variance)
    comp count = count+1
   end

  @PRINCOMP(Ncomp=2) beg tempend+ti-steps pcomp
  # infla1(1) infla1(2) infla1(3) infla1(4) infla1(5) infla1(6) infla1(7) infla1(8) $
    infla1(9) infla1(10) infla1(11) infla1(12) infla1(13) infla1(14) infla1(15) infla1(16)$
    infla1(17) infla1(18) infla1(19) infla1(20) infla1(21) infla1(22) infla1(23)

  IF AVER==0
  {
   set factor1 beg tempend+ti-steps = -pcomp(1)
  }
  ELSE
  {
  set average beg tempend+ti-steps = 0.0

  do i = 1,n
   set average beg tempend+ti-steps = average+(1./n)*infla(i)
  end
  set factor1 beg tempend+ti-steps = average
  }

  do i = 1,n

  if lagopt==1
  {
   do maxlag=1,lag_bic
    linreg(noprint) INFLA_steps(i) beg tempend+ti-steps
    # constant infla(i){steps to steps+maxlag-1} factor1{steps to steps+lag_fact1-1}
    compute bic(maxlag) = log(%rss/%nobs)+log(%nobs)*%nreg/%nobs
   end do

   do maxl = 1,lag_bic
    comp max = %minvalue(bic)
    if max == bic(maxl)
    comp opt_lag(i) = fix(maxl)
   end
  }
  else
  {
   comp opt_lag(i) = lag_fix
  }

  equation eq infla_steps(i)
  # constant infla(i){steps to steps+opt_lag(i)-1} factor1{steps to steps+lag_fact1-1}
  linreg(equation=eq,noprint) INFLA_steps(i) beg tempend+ti-steps
  forecast 1 steps tempend+ti-steps+1
  # eq forint
  set fore_infl(i) tempend+ti tempend+ti = forint
  clear forint

  end do i

 comp inum = inum+1

 END do TI

*** COMPUTE STATISTICS to compare with competitors
*

 do i = 1,n
  DO t = tempend1+steps,end
*   @uForeErrors(noprint) infla_steps(i) fore_infl(i) t-wind t
   @uForeErrors(noprint) infla(i) fore_infl(i) t-wind t
   set rmse_far_serie(i,count_st) t t = %%FERRMSE
  end
 end

************************
** compute statistics RANDOM WALK FORECAST

 smpl beg end

 clear forint

 comp inum = 0

 DO ti = steps,fin

 do i = 1,n
 linreg(define=eq,noprint) INFLA_steps(i) beg tempend+ti-steps
 # infla(i){steps}
 RESTRICT(create,noprint,define=eq)  1 resids1
 # 1
 # 1.0 1
 forecast 1 steps tempend+ti-steps+1
 # eq forint
 set fore_infl_rw1(i) tempend+ti tempend+ti = forint
 clear forint
 end
 comp inum = inum+1

 END do TI


 do i = 1,n
  DO t = tempend1+steps,end
*   @uForeErrors(noprint) infla_steps(i) fore_infl_rw1(i) t-wind t
   @uForeErrors(noprint) infla(i) fore_infl_rw1(i) t-wind t
   set rmse_naiv_serie(i,count_st) t t = %%FERRMSE
  end
 end


**************

** compute statistics SIMPLE AR(p) FORECAST


clear forint
comp inum = 0

DO ti = steps,fin

do i = 1,n

  if lagopt==1
  {
  do maxlag=1,lag_bic
   linreg(noprint) INFLA_steps(i) beg tempend+ti-steps
   # constant infla(i){steps to steps+maxlag-1}
   compute bic(maxlag) = log(%rss/%nobs)+log(%nobs)*%nreg/%nobs

  end do

  do maxl = 1,lag_bic
   comp max = %minvalue(bic)
   if max == bic(maxl)
   comp opt_lag(i) = fix(maxl)
  end
  }
  else
  {
   comp opt_lag(i) = lag_fix
  }
 equation eq infla_steps(i)
 # constant infla(i){steps to steps+opt_lag(i)-1}
 linreg(equation=eq,noprint) infla_steps(i) beg tempend+ti-steps
 forecast 1 steps tempend+ti-steps+1
 # eq forint
 set fore_infl_ar1(i) tempend+ti tempend+ti = forint
 clear forint
 end

 comp inum = inum+1
 END do TI

 do i = 1,n
  DO t = tempend1+steps,end
*   @uForeErrors(noprint) infla_steps(i) fore_infl_ar1(i) t-wind t
   @uForeErrors(noprint) infla(i) fore_infl_ar1(i) t-wind t
   set rmse_ar1_serie(i,count_st) t t = %%FERRMSE
  end
 end

**************

** compute statistics SIMPLE Phillips curve FORECAST
*
clear forint

comp inum = 0

DO ti = steps,fin
*disp ti

do i = 1,n

  if lagopt==1
  {
  do maxlag=1,lag_bic
   linreg(noprint) INFLA_steps(i) beg tempend+ti-steps
   # constant infla(i){steps to steps+maxlag-1} dcp{steps to steps+maxlag-1} ;*ip(i){steps to steps+maxlag-1} m3(i){steps to steps+maxlag-1}
   compute bic(maxlag) = log(%rss/%nobs)+log(%nobs)*%nreg/%nobs

  end do

  do maxl = 1,lag_bic
   comp max = %minvalue(bic)
   if max == bic(maxl)
   comp opt_lag(i) = fix(maxl)
  end
  }
  else
  {
   comp opt_lag(i) = lag_fix
  }
 equation eq infla_steps(i)
 # constant infla(i){steps to steps+opt_lag(i)-1} dcp{steps to steps+opt_lag(i)-1} ;*ip(i){steps to steps+opt_lag(i)-1} m3(i){steps to steps+opt_lag(i)-1}
 linreg(equation=eq,noprint) infla_steps(i) beg tempend+ti-steps
 forecast 1 steps tempend+ti-steps+1
 # eq forint
 set fore_infl_phil(i) tempend+ti tempend+ti = forint
 clear forint
 end
 comp inum = inum+1

END do TI

 do i = 1,n
  DO t = tempend1+steps,end
*   @uForeErrors(noprint) infla_steps(i) fore_infl_phil(i) t-wind t
   @uForeErrors(noprint) infla(i) fore_infl_phil(i) t-wind t
   set rmse_phil_serie(i,count_st) t t = %%FERRMSE
  end
 end

END DOFOR

set far_1 = 0.0
set far_4 = 0.
set far_8 = 0.0
set naiv_1 = 0.0
set naiv_4 = 0.
set naiv_8 = 0.0
set ar_1 = 0.0
set ar_4 = 0.
set ar_8 = 0.0
set phil_1 = 0.0
set phil_4 = 0.
set phil_8 = 0.0

do i = 1,n
set far_1 = far_1+(1./n)*rmse_far_serie(i,1)
set far_4 = far_4+(1./n)*rmse_far_serie(i,2)
set far_8 = far_8+(1./n)*rmse_far_serie(i,3)

set naiv_1 = naiv_1+(1./n)*rmse_naiv_serie(i,1)
set naiv_4 = naiv_4+(1./n)*rmse_naiv_serie(i,2)
set naiv_8 = naiv_8+(1./n)*rmse_naiv_serie(i,3)

set ar_1 = ar_1+(1./n)*rmse_ar1_serie(i,1)
set ar_4 = ar_4+(1./n)*rmse_ar1_serie(i,2)
set ar_8 = ar_8+(1./n)*rmse_ar1_serie(i,3)

set phil_1 = phil_1+(1./n)*rmse_phil_serie(i,1)
set phil_4 = phil_4+(1./n)*rmse_phil_serie(i,2)
set phil_8 = phil_8+(1./n)*rmse_phil_serie(i,3)
end

comp min=1.25, max = 3.
smpl tempend1 end
*spgraph(vfield=1,hfield=3,header='Rolling RMSE: 4 forecast competitors',subheader='average country')
spgraph(vfield=1,hfield=3)
graPH(key=below,min=min,max=max,header='1-step ahead') 4
# far_1
# naiv_1
# ar_1
# phil_1
graPH(key=below,min=min,max=max,header='4-step ahead') 4
# far_4
# naiv_4
# ar_4
# phil_4
graPH(key=below,min=min,max=max,header='8-step ahead') 4
# far_8
# naiv_8
# ar_8
# phil_8
spgraph(done)


/*
open copy rolling_rmse_average.xls
copy(org=obs,for=xls,dates) tempend1 end far_1 naiv_1 ar_1 phil_1 $
far_4 naiv_4 ar_4 phil_4 $
far_8 naiv_8 ar_8 phil_8
close copy
open copy rolling_rmse_far.xls
copy(org=obs,for=xls) tempend1 end rmse_far_serie
close copy
open copy rolling_rmse_ar.xls
copy(org=obs,for=xls) tempend1 end rmse_ar1_serie
close copy
open copy rolling_rmse_rw.xls
copy(org=obs,for=xls) tempend1 end rmse_naiv_serie
close copy
open copy rolling_rmse_phil.xls
copy(org=obs,for=xls) tempend1 end rmse_phil_serie
close copy
*/

**************
/*
DO TA = TEMPEND1+1,END
 DO I = 1,N
  set work_FAR1 I I = rmse_far_serie(I,1)(tA)
 END
 compute frac1=%fractiles(work_FAR1,||0.5||)
 SET FAR_1 TA TA = FRAC1(1)
END

DO TA = TEMPEND1+4,END
 DO I = 1,N
  set work_FAR4 I I = rmse_far_serie(I,2)(tA)
 END
 compute frac4=%fractiles(work_FAR4,||0.5||)
 SET FAR_4 TA TA = FRAC4(1)
END

DO TA = TEMPEND1+8,END
 DO I = 1,N
  set work_FAR8 I I = rmse_far_serie(I,3)(tA)
 END
 compute frac8=%fractiles(work_FAR8,||0.5||)
 SET FAR_8 TA TA = FRAC8(1)
END

***********************

DO TA = TEMPEND1+1,END
 DO I = 1,N
  set work_naiv1 I I = rmse_naiv_serie(I,1)(tA)
 END
 compute frac1=%fractiles(work_naiv1,||0.5||)
 SET naiv_1 TA TA = FRAC1(1)
END

DO TA = TEMPEND1+4,END
 DO I = 1,N
  set work_naiv4 I I = rmse_naiv_serie(I,2)(tA)
 END
 compute frac4=%fractiles(work_naiv4,||0.5||)
 SET naiv_4 TA TA = FRAC4(1)
END

DO TA = TEMPEND1+8,END
 DO I = 1,N
  set work_naiv8 I I = rmse_naiv_serie(I,3)(tA)
 END
 compute frac8=%fractiles(work_naiv8,||0.5||)
 SET naiv_8 TA TA = FRAC8(1)
END

**********************

DO TA = TEMPEND1+1,END
 DO I = 1,N
  set work_ar1 I I = rmse_ar1_serie(I,1)(tA)
 END
 compute frac1=%fractiles(work_ar1,||0.5||)
 SET ar_1 TA TA = FRAC1(1)
END

DO TA = TEMPEND1+4,END
 DO I = 1,N
  set work_ar4 I I = rmse_ar1_serie(I,2)(tA)
 END
 compute frac4=%fractiles(work_ar4,||0.5||)
 SET ar_4 TA TA = FRAC4(1)
END

DO TA = TEMPEND1+8,END
 DO I = 1,N
  set work_ar8 I I = rmse_ar1_serie(I,3)(tA)
 END
 compute frac8=%fractiles(work_ar8,||0.5||)
 SET ar_8 TA TA = FRAC8(1)
END
***********************

DO TA = TEMPEND1+1,END
 DO I = 1,N
  set work_phil1 I I = rmse_phil_serie(I,1)(tA)
 END
 compute frac1=%fractiles(work_phil1,||0.5||)
 SET phil_1 TA TA = FRAC1(1)
END

DO TA = TEMPEND1+4,END
 DO I = 1,N
  set work_phil4 I I = rmse_phil_serie(I,2)(tA)
 END
 compute frac4=%fractiles(work_phil4,||0.5||)
 SET phil_4 TA TA = FRAC4(1)
END

DO TA = TEMPEND1+8,END
 DO I = 1,N
  set work_phil8 I I = rmse_phil_serie(I,3)(tA)
 END
 compute frac8=%fractiles(work_phil8,||0.5||)
 SET phil_8 TA TA = FRAC8(1)
END


smpl 2000:1 end
graph(key=below) 2
# infla_steps(2)
# fore_infl(2)
# fore_infl_rw1(2)
# fore_infl_ar1(2)
# fore_infl_phil(2)
*/









